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Light bosonic degrees of freedom have become a serious candidate for dark matter, which seems 
to pervade our entire universe. The evolution of these fields around curved spacetimes is poorly 
understood but is expected to display interesting effects. In particular, the interaction of light 
bosonic fields with supermassive black holes, key players in most galaxies, could provide colourful 
examples of superradiance and nonlinear bosenova-like collapse. In turn, the observation of spinning 
black holes is expected to impose stringent bounds on the mass of putative massive bosonic fields 
in our universe. 

Our purpose here is to present a comprehensive study of the evolution of linearized massive scalar 
and vector fields in the vicinities of rotating black holes. The evolution of generic initial data has a 
very rich structure, depending on the mass of the field and of the black hole. Quasi-normal ringdown 
or exponential decay followed by a power-law tail at very late times is a generic feature of massless 
fields at intermediate times. Massive fields generically show a transition to power-law tails early on. 
For a certain boson field mass range, the field can become trapped in a potential barrier outside the 
horizon and transition to a bound state. Because there are a number of such quasi-bound states, 
the generic outcome is an amplitude modulated sinusoidal, or beating, signal, whose envelope is well 
described by the two lowest overtones. We believe that the appearance of such beatings has gone 
unnoticed in the past, and in fact mistaken for exponential growth. The amplitude modulation of 
the signal depends strongly on the relative excitation of the overtones, which in turn is strongly tied 
to the bound-state geography. 

A fine tuning of the initial data allows one to see the evolution of the nearly pure bound state 
mode which turns unstable for sufficiently large black hole rotation. For the first time we explore 
massive vector fields in generic BH background which are hard, if not impossible, to separate in 
the Kerr background. Our results show that spinning BHs are generically strongly unstable against 
massive vector fields. 

PACS numbers: 98.80.Es,11.25.Wx,14.80.Va,04.70.-s 
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I. INTRODUCTION 

One of the most exciting outcomes of General Relativity (GR) are black holes (BHs), the physics of which has 
grown into a mature and fully developed branch of GR and extensions thereof PQ[2]. Observations of, e.g., X-ray 
binaries indicate that solar mass (3 — 30Mq) BHs mark the endpoint of the life of massive stars and are anticipated 
to be a significant component of the galaxies' population. Supermassive BHs (SMBHs) with masses 10^ — 10^ Mq or 
higher are conjectured to be hosted in the center of most galaxies, controlling galaxy growth and evolution, stellar 
birth and powering active galactic nuclei and other powerful phenomena. 

Tremendous progress has been made in actually observing some of the fascinating general relativistic effects. From 
X-ray spectra on the inner edge of accretion disks, which probe the innermost stable circular orbit of the geometry, 
to gravitational wave (GW) physics, "precision BH physics" is a new and rapidly developing field |^{5]. The future 
holds the promise to observe some of these effects accurately by monitoring the supermassive BH at the center of our 
own galaxy. 

One of the fundamental reasons why precision BH physics is possible at all, are the no-hair and uniqueness theorems: 
BHs in 4-dimensional, asymptotically flat spacetimes must belong to the Kerr-Newman family and are, thus, fully 
specified by three parameters only: their mass, angular momentum and electric charge (see e.g. Ref. |6','7j, or Carter's 
contribution to Ref. 0). In more colloquial terms, this is commonly expressed by saying that BHs have no hair or, 
rather, have three hairs only. This simple yet powerful result has far reaching consequences: Given some arbitrary 
perturbations with the same conserved charges, they must all decay to the same final state, namely one BH with those 
charges. By now, there are a plethora of studies, at the perturbative and fully non-linear level, investigating how this 
unique final state is approached (see, e.g., Ref. pi 110) for recent overviews). In the following we briefly summarize 
these studies. 
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FIG. 1. Time evolution of a dipole (l — l,m = 0) scalar Gaussian wave packet in Schwarzschild background. We clearly 
observe the main features of such a field: (i) a prompt response at early times followed by (ii) the quasi-normal mode ringdown 
and (iii) a late-time tail. 

Generic response of a BH spacetime to external perturbations. The generic behavior of massless fields 
around a BH is illustrated in Fig.[lJ where we plot the evolution of a Gaussian wave packet \1/ = e"^'"^"-* around 
a Schwarzschild BH. The particular initial data refers to a scalar field, but the qualitative results are universal and 
independent of the initial conditions. The generic behavior of massless fields around a BH can be divided into three 
parts (c.f. Fig. [T]): 

(i) An initial data-dependent prompt response at early times, which is the counterpart to light-cone propagation in 
flat space; 

(ii) An exponentially decaying "ringdown" phase at intermediate times, where the BH is ringing with its characteristic 
quasi-normal modes (QNMs). This stage typically dominates the signal, and its properties, such as vibration frequency 
and decay timescale, depend solely on the parameters of the final BH '3] . Because of the no-hair theorem, the detection 
of QNMs allows one to uniquely determine the BH mass and spin and provides tests of GR ^ [TTJ [T^ ; 

(iii) At late times, the signal is dominated by a power-law fall-off, known as "late-time tail" [OVllSj . Tails are caused 
by backscattering off spacetime curvature and more generically by a failure of Huygen's principle. As such, tails also 
appear in other situations where light propagation is not on the light cone such as in massive field propagation in 
Minkowski spacetime |16j . or massless field propagation in odd-dimensional spacetimes |17j . 

Superradiant effects. The long-lasting oscillation of the lowest QNMs is the most important stage in the life of any 
field around a BH. Its lifetime, or quality factor, depends solely on the BH spin p]. Specifically, the lifetime tends 
to increase with growing spin and the decay timescale approaches zero for nearly extremal BHs. This behaviour is 
tightly connected to superradiance |18H20) : In a scattering experiment of low- frequency waves off a BH the scattered 
wave is amplified if the real part of its frequency ujji satisfies the superradiant condition 

ujR < mVtH , (1) 

where m is the azimuthal "quantum" number and U,h is the angular velocity of the BH horizon. We refer the reader 
to App. |A] for a derivation of this condition for both, scalar and vector fields. The excess energy is withdrawn from 
the object's rotational energy [TH1II2] and, in a dynamical scenario, the BH would presumably spin down. The effect 
can be attributed to the existence of negative-energy states in the ergo-region, and dissipation at the event horizon. 
Superradiance is the chief cause of a number of exciting phenomena in BH physics: 

(i) Generic perturbations are damped away to infinity and across the event horizon. Because rotating BHs amplify 
waves that fulfill the superradiant condition, Eq. ([T]), the amplification factors as well as the quality factor of these 
superradiant modes increase with rotation. 

(ii) Satellites around BHs typically spiral inwards as time goes by, due to gravitational wave emission and energy 
conservation. Emission of radiation to infinity results in a larger binding energy of the particle. Because superradiance 
implies the extraction of the BH's rotational energy, it is possible that the energy deficit comes entirely from the BH 
kinetic energy. In this way, satellites around rapidly spinning BHs can in principle orbit at a nearly fixed radius - 
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on so-called floating orbits - for a much longer time, tapping the BH's kinetic energy. In BH binaries, this effect can 
dominate in the presence of resonances pm23j . This phenomenon is analogous to tidal acceleration, e.g., in case of 
the Earth- moon system [531 US] . 

(iii) A further interesting effect can be triggered by enclosing the spinning BH inside a perfectly reflecting cavity. As 
was recognized already by Zel'dovich (18i il9j , any initial perturbation will get successively amplified near the BH and 
reflected back at the mirror, thus creating an instability, which was termed the "BH bomb" [261 . Whereas the 
setup appears physically artificial at first glance, the role of the mirror can actually be realized naturally in many 
ways, including an anti-de Sitter spatial infinity. In this case, the BH bomb translates into a real, physical instability 
of (small) rotating BHs in asymptotically AdS spacetimes [5SH3T]. 

(iv) Finally, of direct interest for the present study is the fact that massive fields around Kerr BHs are also prone to 
a BH bomb-like instability, because the mass term effectively confines the field [551 - 155] . 

Consider a scalar field surrounding a black hole with mass M and angular momentum J = aM . The instability is 
regulated by the dimensionless parameter fisM (from now on we set G — c = 1), where TOj = fisfi is the scalar field 
mass, and is described by the time dependence of the field, ^ ^ g-i^t -^j^j-^ complex frequency uj = uir + lOJi. For 
small coupling M^s <^ 1 the characteristic (unstable mode) frequency giving rise to the instabihty is [iOlliT] 

In the opposite limit, i.e., for very large mass couplings M/ig >> 1, the characteristic inverse time is [33] 

Mw/ =10-^exp(-1.84Af/is) ■ (3) 

The instability timescales are typically large. The scalar field growth rate has a global maximum of r = l/w/ ~ lO^M 
for the dipole with mass coupling Mfjis = 0.42 in the background of a Kerr BH with a/M = 0.99 [36] l37]. 

The above results refer to massive scalar fields in the background of Kerr BHs. It was widely believed that massive 
vectors would be subject to a similar instability. Unfortunately, the non-separability of the field equations renders 
this is a non-trivial problem. Recently, significant progress has been made, with a thorough study of massive vector 
fields around Schwarzschild BHs and slowly rotating Kerr BH backgrounds [lOlllI]- Pani et alnse a slow-rotation 
expansion of the Fourier transformed field equations, accurate to second order in rotation, to prove that the Kerr 
spacetime is indeed unstable against massive vector fields [IHl HI] • The massive vector field instability can be orders 
of magnitude stronger (i.e., shorter timescales) than its scalar counterpart. 

All these calculations have been performed in the linear regime, thus neglecting backreaction effects such as the BH 
spin-down or effects due to non-linear self-interaction of the scalar field. Therefore, the final state of the superradiant 
instability in the fully dynamical regime is not known, partly because it requires the non-linear evolution of Einstein's 
equations for a timescale of order lO^M. A plausible evolution scenario consists of an exponentially growing scalar 
condensate outside the BH, extracting energy and angular momentum from the BH until the superradiant extraction 
stops, i.e., until the condition ([T]) is no longer satisfied. Further interesting new phenomena arise when we consider 
non- linear interaction terms, such as bosenova-type collapse presented in Refs. [35H3S]; or higher dimensional back- 
ground spacetimes, such as the boosted black string recently reported in [46] or the Schwarzschild-Tangherlini solution 
discussed in [ITlliB]. 

Superradiant instability in astrophysical systems. Massive fields in the vicinity of BHs are subject to a 
BH bomb-like, superradiant instability, and they grow exponentially with time. However, the effect is very weak for 
known standard model particles in astrophysical environments: For example, the mass coupling for the lightest known 
elementary scalar particle, the pion, around a solar mass BH is M^s ~ 10^®, resulting in an instability timescale 
much larger than the age of the universe. Nevertheless, the superradiant instability might become significant if we 
consider standard model particles around primordial BHs (see, e.g., |49n50n or if there exist fields with small, but non- 
vanishing mass. One exciting possibility for these fields is provided by axions, ultralight bosonic states emerging from 
string-theory compactifications, which have not been ruled out by current experiments. In the "axiverse" scenario 
an entire landscape of ultra-light pseudo-scalar fields covering a mass range from 10~^^eV < /is < lO^^eV^ has 
been proposed (see [51 [4] [43] for recent overviews). In fact, the existence of ultra-light axions leads to a plethora 
of possible observational implications and signatures, such as modifications of the cosmic-microwave background 
polarization (for 10^^^ eV < /Lts < 10~^^eV). They are also anticipated to make up a fraction of dark matter if 
10~'^^eV < /xs < 10~^^eV 5", "WT, "5? . Of particular interest in the context of BH physics are axions in the mass range 
10~^^eV < /xs < 10~^°eF [3, 4, 43 . Then, the time scales for the superradiant instability become astrophysically 
significant, giving rise to a number of interesting effects: 



^ Notice the difference of a factor 2 to the original result | 34 |. 
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(i) A bosonic cloud bounded in the vicinity of a Kerr BH might create a "gravitational atom" , which can be de-excited 
by the emission of gravitons, thus carrying away BH angular momentum; 

(ii) If the accretion of bosons from this cloud is efficient enough, the rotation of the BH can be sustained and it might 
be turned into a GW pulsar; 

(iii) If, on the other hand, the accretion from the axionic cloud is not efficient enough, the BH will eventually spin 
down, thus yielding gaps in the Regge plane (the phase-space spanned by mass and spin parameter of the BH). Further 
possible effects have been discussed in Refs. HH HI IHSSl [S3]. 

Similar superradiant instabilities are expected to occur for massive hidden [/(I) vector fields, which are also a 
generic feature of extensions of the standard model [MH57] . As already stated, while superradiant instabilities have 
been widely studied for massive scalar fields [211 EH EH ESI |36l |37l EHl [59] , the case of massive vector fields is still in 
its infancy, though significant work along these lines was recently reported flOHi^ ITTl [15] . 

So far most studies on the massive boson instability have been performed in Fourier space. An early attempt 
at studying the massive scalar field instability in the time domain, with generic initial conditions was presented by 
Strafuss and Khanna [6^. We believe that, while the technical study may be correct, some of its conclusions are not; 
specifically, the authors reported an instability growth rate of Mojj ~ 2- 10~^, which is two orders of magnitude larger 
than previous results in the frequency domain ^36i i37j and more recent numerical studies in the time domain [44i, i6l] . 
We will attempt a correct explanation for these puzzling results in the body of this work. 

The purpose of the present study is to investigate the time evolution of generic linearized massive scalar and vector 
fields in the vicinities of spinning BHs. Surprisingly, not much seems to have been done on this problem. Our 
"generic" initial data consists of Gaussian wave packets, but we will also study the evolution of bound state modes. 
The exploration of nonlinear gravitational dynamics or self-interactions will be presented elsewhere. 

This work is organized as follows: In Sec. [H] we present the numerical framework, describing the formulation as 
a Cauchy problem, the setup of initial configurations and the background spacetime. Sec. |III| is devoted to the 
numerical results of massive scalar field evolution. In particular, we present a number of benchmark tests to verify 
our implementation before studying more generic setups. We will show that the evolution of a massive scalar has 
a non-trivial pattern, which can be explained in terms of multi-mode excitation. We believe that this pattern also 



describes the results reported by Strafuss and Khanna [60^. In Sec. IV we discuss our investigations of the massive 
vector (also known as Proca [57]) field in generic Kerr BH backgrounds, where we show, for the first time in a time 
evolution of rapidly spinning BHs, that Kerr BHs are strongly unstable against these fields. Finally, we summarize 
our results and present concluding remarks in Sec. jVj 

II. SETUP: ACTION, EQUATIONS OF MOTION AND BACKGROUND METRIC 

A. Action and equations of motion 

We consider a generic action [S] |62l involving one complex, massive scalar ^ and a massive vector field with 
mass ms — fJ-sfi and niv — fJ-vfi, respectively. 

Here, the potential V'(vE') is of cubic or higher order in the scalar field. The scalar and vector fields are allowed to 
interact through the axion-like coupling constant fcaxion- Ffii^ = V^Ajy — V,yA^ is the Maxwell tensor and * F^'^ = 
i^t^iypap ^^^i jjgj.g^ gMi'p'T ^ _iptj.i^pa j^j^^ ppi^P'y is the totally anti-symmetric Levi-Civita symbol with 

^0123 — 2_ |2j ^j^g resulting equations of motion are 

(V.V^ - 4) * = ^ ^F^^'F^. + V'i^) , (5a) 
^.Ff^" + My^^ = -2fcaxio„ *F^''9.* , (5b) 

fi* a*'" + + V{^) \ + \ {^*■^'^''' + 4"^^'*.'') . (5c) 




2-^ V2 2 ' ' A 



The identity Vi^ *F^^" = is useful to derive the equations of motion for the Chern-Simons term. 
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We note that these equations describe the fully non-linear evolution of the system. Also, we have written the equations 
such that all terms quadratic or of higher order in the vector or scalar fields appear on the right hand side. In the 
remainder of this work, we will restrict ourselves to the case of scalar and vector fields with small amplitudes and will 



ignore the higher-order contributions on the right-hand sides of (5a I- (5c 



Under this assumption Eq. (5c) is equivalent to Einstein's equations in vacuum and a solution to this equation is 



the Kerr metric which in Boyer-Lindquist coordinates is given by 



ds^ = - 



2Mr 



BL 



AMr 



BL 



dtdr 



BL 



AMrsLa sin^ 9 



^d0^ + sin^ 
dtd(l) - 2a sin^ 9 ( 1 



BL + a 
2A/rBL 



2Ma^rBL sin^ 



d(j)^ 



dr 



BLC 



with 



' BL 



cos^ 9, A 



'BL 



2Mr 



BL 



(6) 



(7) 



This geometry describes a rotating BH with mass M and angular momentum J — aM . Note that in order to ensure 
the regularity of the spacetime, i.e. the existence of an event horizon, the BH spin is constrained by the Kerr bound 
a/M < 1. 

A second consequence of our assumptions is that the axionic coupling can be neglected. This means that we 
effectively study minimally coupled massive scalar and vector fields separately and our results will describe small 
linearized fields around the Kerr background. Any potential instability we find is consistent with the above assumptions 
for timescales small enough such that the fields are small. Over long timescales, the fields may grow to large amplitudes 
where our assumption no longer remains valid and a non-linear study becomes necessary. We postpone such non-linear 
evolutions to a future investigation. 



In our approximation, the scalar and vector field dynamics are governed by the linearized version of Eqs. (5a 
and (5b) 



(V.V-mDvI'^O, 



(8a) 
(8b) 



while the Kerr metric ^ satisfies G^^, = 0, i. e. Eq. (5c 



linearized in and A, 



Evolution equations for scalar fields. Because we intend to solve the equations of motion iSan and (8b) 



numerically, it is convenient to reformulate them as time evolution problem. For this purpose we employ the 3 + 1- 
decomposition of the spacetime (see e.g. [SSj) and consider the background spacetime in generic 3 -|- 1-form 



ds^ = - a^dt^ + Jrjidx' + P'dt){dx^ + P^dt) 



(9) 



Here, 7ij is the spatial metric and a and /3* are the lapse function and shift vector which represent the coordinate or 
gauge freedom of general relativity. We introduce the conjugated momenta 



ni = --{dt 
a 



(10) 



where Xji := 3?(X) and Xj := '^(X) denote the real and imaginary parts, respectively. Definition (10) provides 
evolution equations for the scalar field 4" 



aHf 



dt'^i = Cp'^i - aH/ , 



where Cf^'^ rj — P^dk"^ rj. By applying the 3 + 1-split to the Klein-Gordon equation (8a 
equations for the momentum 

9tHfl ^CfiWR - D'aD.^R + a{-DW,^R + KHr + 4*«) 



(11) 

we obtain the evolution 

(12a) 
(12b) 



where Ci^IIrj = p^diJlRj. Di is the covariant derivative associated with the 3-metric 7^ and K is the trace of the 
extrinsic curvature. 



7 



Evolution equations for vector fields. We next apply the 3+1 decomposition to the evolution equation (8b) and 
obtain 

VV^A, - VV.A,, + filrA,, = - [V^V.A^ - V^IVA,) - Rf.^'A, - fi^A^] = . (13) 



By operating with on Eq. ( 13 ) it is straight-forward to show that the Lorenz gauge 

V^A^ =0 (14) 
needs to be satisfied. For a vacuum background spacetime as considered in our work, we also have R^i, = and 



Eq. ( 13 ) simplifies to 

VV.A,, - fi^A,, = . (15) 



Note, that in case of a non- vanishing cosmological constant A, i?^^ — A^^^ and in place of Eq. ( 15 ) we would obtain 



V'V,A^-iA + ^ii,)A^^O, (16) 

i. e. the cosmological constant enters as an additional "mass" -like term. In particular, it changes the evolution equation 
of a massless vector field to that of a massive one. Note, however, that the Maxwell equations in Kerr-(anti-)de Sitter 
background are known to be separable and can be written in the form of a Teukolsky type equation [2^1 ISH HI] ■ 
Therefore, one might expect that the equations of motion for a massive vector field in a vacuum Kerr spacetime 
should also be separable. We emphasize, however, that this analogy between a massless field with cosmological 
constant and a massive one without needs to be taken with care: a massless vector field has only 2 dynamical degrees 
of freedom irrespective of the value of the cosmological constant, whereas a massive vector field has 3. In fact, up 
to date, a separation of the equations of motion for a massive vector field has not been accomplished. It is not 
immediately obvious, for this reason, whether there exists a well-defined correspondence between the two cases. In 
this work we focus on A = and therefore leave a detailed investigation of this question for future work. 

We now apply the 3-1-1 decomposition to the vector field and split A^ into its spatial part and normal component 

A' = t'^m A ' and (p^ -ni'A^, , (17) 

where is the vector normal to the spatial hypersurface S. The vector field can be reconstructed from its projections 
according to A^ = -\- n^ip. Furthermore, the projection of the Maxwell tensor along the normal vector yields 
the electric field 

Ef, ^F,,,rf , (18) 

which is a purely spatial quantity, i.e., E^n^ = 0. 

With all the necessary ingredients at hand we now proceed by performing the 3 + 1-split Eqs. (14) and (15). In 
terms of the dynamical variables {ip, At, Ei} this procedure results in the constraint 

Ce ^D'E, + //^(^ = . (19) 

and in the evolution equations 

[dt - Cf})ip = - a {D'A^ - Kip) - A^D'a , (20a) 
[dt - Cp)A^ ^~a{E, + D,p>) - p>D,a , (20b) 
{dt ~ Cp)E, =a {^j.lA, + KE, - 2E^K,j + {D.A, - D^A,)) + D^a{D,A, - D^A,) , (2Gc) 

where Cpp> = P^dkP>, CpA^ = P^dkA^ + Ak^^P^ and CpE, = p^dkE, + Ekd.p^ . 



B. Background in horizon penetrating coordinates 

In practice, it is convenient to employ horizon penetrating coordinates and consider the Kerr spacetime in Cartesian 
Kerr-Schild coordinates {t,x,y,z). Without loss of generality, we assume the angular momentum to point in the z 
direction. Then, the Kerr-Schild form of the lapse function a, shift vector , 3-metric and extrinsic curvature 
is given by 

a={l + 2HlH')-''\ p^ = -——^^, j,,=5,,+2Hkl,, (21a) 
K,, = - -ihljdtH + 2Hl(,dtlj)) - 2a (%(?j)i/;*) + 2HHU%,d\k\lj) + HlH,ljl''dkH) , (21b) 
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where 

Mrl^ _ f tblX + ay tblV - ax z \ 

and the Boyer-Lindquist radial coordinate tbl is related to the Cartesian Kerr-Schild coordinates by 

^' + ^^' + — =1 (23) 

C. Initial data 

In this work we consider two types of initial configurations: (i) generic pulses of Gaussian shape and (ii) bound 
states which are particularly suitable for identifying putative instabilities. We describe each of these initial data in 
detail. 

Gaussian initial data. We specify Gaussian wave packets of the form 

*(t = 0)=0, n(t = 0)=cxp(^-^'^^^^ oS(0», (24a) 

(^(t = 0)=0, A(i = 0)=cxp(^-^'l— _iS(0,0), E,it = 0) = 0, 2 = 1,2,3, (24b) 

where r = -^/x^ + + is the Kerr-Schild radial coordinate, tq and w are the center and width of the Gaussian, 
while oT,(6,(f)) and _il](0,0) represent superpositions of spherical harmonics sYim{Q,4') of spin weight s = and 
s = —1, respectively. Expressed in Cartesian coordinates [x^y^z) 

X — r sin 6 cos (/) , y — r sin 6 sin , z = r cos 6 , (25) 



the spin- weighted spherical harmonics up to Z = 2 are given by Eqs. ( |B1[ )-(B6 ) in Appendix |B[ 

Bound state initial data. Our second type of initial data is given by the quasi-bound states of massive scalar fields 
around BHs which represent long-lived modes of massive scalar field perturbations around Schwarzschild or Kerr BHs 
and have been studied extensively in the literature in the frequency domain [SJ [Ml [351 133 HOI HIl US] ■ Our particular 
interest in these modes arises from their pure nature; they represent potentially superradiant, single- frequency states. 
By specifying single-mode states of this type, we are able to suppress interference or beating effects of the kind 
discussed below for Gaussian initial data. Evolutions of such modes additionally serve as a useful test for our code 

m- 

There exist powerful and simple methods to construct the bound states for massive scalars, either by direct numerical 
integration or a continued fraction approach [9l [HI |36l [37] . Here, we adopt Leaver's continued fraction method and 
obtain in Boyer-Lindquist coordinates 

^irn =e-^"*-^e-™'^--5z™(0BL)i?im(rBL) • (26) 
Here, ^'/^(^bl) ^i'G spheroidal harmonics [67 and the radial dependence is given by 

OO ^ -.71 

RUtbl) =(rBL - rBL,+)-^'^(rBL - rBU-r+^-'e^^-" V a„ ( ""^^ ~ '"^^'+ ) , (27) 

with 

a= , q = ±Jni;-uj^, x=— ■ (28) 

rBL,+ -rBL,- V q 

rBL,± — M ± V Af 2 — are the radii of the inner and outer horizon and uic = raVLn = fnjMr^ — ^^'^ critical 
frequency for superradiance. All remaining terms in this expression are known in closed form and the characteristic 
frequency ui can be obtained by solving a three-term recurrence relation for the coefHcients a„ given by, e.g., Eqs. (35)- 
(48) of [37]. For our purposes, we still need to transform these results from Boyer-Lindquist to Kerr-Schild coordinates 
denoted in this discussion for clarity by a subscript "KS" , 

dtKS ^dtBL + ^^^BL ^^^^ ^ ^^^^ ^ ^^^^ ^ ^^^^ ^ ^^^^ ^ ^^^^ ^ ^^^^ ^ ^drBL • (29) 
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Then, the bound state scalar field is given by 

*;m =e "^^(r-KS - '^KS.+j (»'ks - ^-ks,-) 

V^KS-?'KS-/ 



^;m(^KSi 0Ks)-R/m , (30) 



where A = — 2'"*'^'-ks,+ ^ ^ 2»a;MrKs,- ^ . Henceforth we will drop the subscript "KS" and denote 

»'KS,--»'KS,+ ' '•KS,--'"KS,+ ' '•KS,--''KS,+ ^ ^ 

the Kerr-Schild radius and angles by (r, 9, (j)). 



The corresponding conjugated momenta Him are computed from their definition, Eq. (10 1. Unfortunately, a simple 
construction of these modes exists only for scalar fields, while for vector fields a fully numerical procedure is required 
[551 1571 HOI HT] , Partly for this reason and partly because Gaussian initial data turn out to be adequate for the 
identification of superradiant instabilities we will not consider vector bound states in the remainder of this work. 



D. Wave extraction and output 



The main diagnostic quantities extracted from our simulations are the radiated scalar and vector waves. The scalar 
multipoles are directly obtained from interpolating the fields and 11 onto spheres of constant coordinate radius 
r — and projecting onto 5 = spherical harmonics according to 



J dn^iJit, d, 4^)Yi;^{9, ^) , nut) = J dnnit, 0)ii:„(0, 0) . (31) 



For vector fields, we construct the gauge invariant Newman Penrose scalar [BS] (see also, e.g., [5M7T] for recent 
applications in numerical simulations). 

*2 ^F^J'-fh'' , (32) 

where — -^{n'^ — u^) and = ^(^^ — tw^) are vectors of a null tetrad. In practice, the vectors of the null 

tetrad are constructed from a Cartesian orthonormal basis on the spatial hypersurface and the timelike 

orthonormal vector n'^. Together with the reconstruction of the Maxwell tensor from 

Ff_,i, ^Hf^Ei, - n^Efj_ + Df^Au - Di,A^ (33) 

we can straightforwardly derive the Newman-Penrose scalar $2 whose real and imaginary components are given by 

= - ^ [E^v' + u'v' {D,Af - DjAf) + El w' + u'w' {D,A] - DjAj )] , (34) 
*2 = I [eS' + (D^Af - D,A^) - El - u^v^ {D,Al - D,AD] . (35) 
Finally, we obtain the multipoles by projecting $2 onto s = — 1 spin- weighted spherical harmonics 

KJi)^ j dn[^^{t,eA)-iYil{o,<l^) + ^i{t,e,^)^,Yi,{e,^)\ , (36) 



E. Numerical implementation 

We have implemented this framework as a module Lin-Lean in the Lean code [75], which is based on the Cac- 
tus computational toolkit [751 [H] and the Carpet mesh refinement package [7S1 [7S]. The evolution equations 
are integrated in time using the method of lines with a fourth-order Runge-Kutta scheme and fourth-order spatial 
discretization on all refinement levels except the innermost where we excise the black-hole singularity and use second- 
order accurate stencils instead. The excision is implemented by removing a "legosphere" of radius r^^ < IM centered 
on the singularity [77| [78] and therefore guaranteed to be inside the event horizon for all values of a. Inside the 
excision region, we set all evolution variables to their flat spacetime values a = 1, = 0, ^ij — Sij and Kij = 0, 
and extrapolate values onto the excision boundary from the exterior regular grid; see [77[ 178] for more details of this 
procedure. 
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Run 




a/M 








Grid Setup 


sS_mOOO 


0.00 


0.00 


Voo 


+ Yio + yi-i 




{(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/lOO} 


sK_mOOO 


0.00 


0.99 


YoQ 


+ yio + yi-i 


- Yii 


{(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = A//108} 


sS_m001 
sS_m042 
sS_mlOO 


0.10 
0.42 
1.00 


0.00 
0.00 
0.00 


Yio 


Vii 

+ Yii + Y20 + Y22 
Yii 


{(1024,512,256,128,64,32,8,4,2), h = M/40} 
{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 
{(1024, 512, 256, 128, 64, 32, 8, 4, 2), 1/40} 


sK_m042c 

sK_m042,„ 

sK_m042/ 


0.42 

0.42 
0.42 


0.99 
0.99 
0.99 


Yoo 


+ Yio + Yi-i 
Vii.i-i 
5^11,1-1 


-Yii 


{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 
{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/72} 
{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/84} 



TABLE I. Initial setup for simulations of a scalar field with Gaussian initial data located at ro — 12 M and with width 
w = 2M in a Schwarzschild or Kerr background. We denote the mass parameter M^s, the dimensionless spin parameter a/M, 
the modes of the initial pulse qT,{9, (j)) and the specific grid setup, measured in units of the BH mass M , following the notation 
of Sec. II E in [72|. 



M/.IS 


a/M Mcjii (n 


= 0) 


Mojii (n 


= 1) 


M5w 


0.00 


0.00 


0.2929 - 


i0.09766 


0.2645 - 


i0.3063 


0.0284 


0.00 


0.99 


0.4934 - 


i0.03671 


0.4837 - 


i0.09880 


0.0097 


0.42 


0.00 


0.4075 - 


iO. 001026 


0.4147- 


i0.0004053 


0.0072 


0.42 


0.99 


0.4088 + 


il.504- lO"'^ 


0.4151 + 


i5.364 ■ 10"* 


' 0.0063 


1.00 


0.00 


0.9222 - 


iO. 11842 


0.9570 - 


i0.04792 


0.0348 


1.00 


0.99 


0.8765 - 


iO. 18994 


0.9515- 


i0.06292 


0.075 



TABLE II. Frequencies for the fundamental (n = 0) and first overtone (n = 1) modes of (i) the quasi-normal modes for 
M^s ~ and (ii) the bound states for Alfis = 0.42, 1.0 for a set of BH spins a/M obtained with the continued fraction 
method. The beating frequency Sio is defined as the difference of the real parts of the frequencies for n — and n — 1. 

III. SCALAR FIELD EVOLUTIONS 

In this section we report our results obtained for the evolution of scalar fields in BH background spacetimes. For 
this purpose, we have evolved Gaussian wave pulses of width w = 2 M centered at ro = 12 M around either a non- 
rotating Schwarzschild BH or a rapidly spinning Kerr BH with a/M = 0.99. The set of our simulations is summarized 
in Table U 

As mentioned in the introduction, the resulting signal is typically composed of three stages, a transient, an expo- 
nential ringdown phase and late-time tails. Whereas we identify these stages clearly in evolutions of massless scalar 
fields, the massive case exhibits a richer phenomenology which we will discuss in detail further below. First, however, 
we summarize results from the literature and present tests of our numerical infrastructure. 

A. A summary of results in the literature 

In the following discussion, we will make frequent use of the characteristic frequencies of the two lowest (fundamental 
and first overtone) modes of the dipole as obtained from a linearized analysis in the frequency domain [21 |TT] . These 
values are summarized for a set of chosen mass parameters /is of the scalar field and our two choices of the rotation 
parameter a/M in Table [lT| 

In the case of massless perturbations the classification of oscillation modes is relatively straightforward; there 
exists one family of oscillation modes, the quasi-normal modes, characterized by an integer number and the "lowest" 
or fundamental mode is defined as having the largest damping time. In fact, massless perturbations are generally 
short-lived and decay fast (on timescales of order of the BH mass) unless the BH is rotating at nearly extremal 
rate. In contrast, massive perturbations generally decay on larger timescales and their class of oscillation modes 
contains a second family with profiles concentrated near the effective potential well. These long-lived modes, referred 
to as "bound states" |37j . are conventionally ordered by decreasing absolute value of the imaginary part, i. e. the 
fundamental oscillation mode is defined as that with the shortest damping time. This is in contrast to the usual 
quasinormal modes, which are conventionally ordered by increasing imaginary part. The values for the fundamental 
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(n = 0) and the first overtone (n = 1) bound state mode obtained from the continued fraction method are hsted in 
Table |ll| together with the difference Siq = oj^'^^ —uj^^ which wiU play an important role in our interpretation of the 
resultsiurther below. 

In the regime of small mass couplings Mfj,s, one finds an analytic approximation for the oscillation frequencies of 
these bound states, resembling a hydrogen-like spectrum [iOl l¥T] . 



2{l + n + iy 



M.,^-...(^)^ (38b) 

where Siy is given by Eq. (C8) of [41]. 

At late times, after the relaxation of the BH, power-law tails arise [T3HI51[Zn]- Power-law tails exist both in flat 
and in curved spacetime and are generically caused by a failure of Huygen's principle. In the case of flat spacetime, 
they can arise in massive interactions [16] or even with massless fields in odd-dimensional spacetimes due to the 
peculiar nature of the Green's function 117'. Curved backgrounds give generically rise to power-law tails which are a 
consequence of the continued scattering of the scalar field off the background curvature. These tails have the form 

* - sm{fj.st) . (39) 

For massless scalar fields, one finds one exponent p given by 

p=-{2l + 3). (40) 

Tails of massive perturbations, however, exhibit a more complex behaviour; details depend on the mass parameter fis 
[80H83j ■ but their intermediate and late-time behaviour is characterized by 

p = — {I + 3/2) , at intermediate times , (41a) 
p — —5/6 , at very late times . (41b) 

These exponents are loosely associated with the two timescales M, fj,s of the problem; in particular, intermediate 
times refer to the window (we assume small M ns) M^s <C t^s ^ (M/xg)"^. 

Note, that the intermediate-time behavior is identical to the power-law behavior of massive fields in flat spacetime 
[l6] whereas the late-time behavior is also affected by scattering off the spacetime curvature. 

B. Space dependent mass coupling and massless scalars 

We next employ the results summarized in the previous section to test our numerical framework. For the first 
test, we consider the unphysical scenario of a scalar field with space-dependent mass term /i| = — lOM^/r'* in a 
Schwarzschild background. This choice is motivated by the strong instability it generates and provides a unique and 
fast setup to test the code in a particularly violent regime. A mode analysis of the Klein-Gordon equation for this 
choice of us is straightforward and demonstrates the existence of at least one unstable mode with time dependence 
xj; gO. 071565*^ 'pjjg results obtained for the evolution of a spherical shell described by '^t=o — ^/l^ ^^P ^'^"l^"' ^ 
are shown in Fig. [2] which shows the amplitude of the scalar field extracted at rex = 10 M as a function of time. Our 
numerical results are consistent with an exponential growth, ^' ^ g0 07i6it excellent agreement with the frequency 
calculation. 

A second test for our code is provided by the evolution of massless scalar fields around BHs, a case well studied 
in the literature [SKTUI. For this purpose we have initialized the field by a Gaussian with rg = 12 M and w = 2 M 
and extracted the monopole and dipole at r^x = 10. The results are shown in Fig. |3] for a BH background with 
a/M — and a/AI ~ 0.99, respectively. The waveform displays the familiar features; an early transient followed by 
an exponentially decaying sinusoid and a power-law tail at late times. A fit to the ringdown phase of the dipole yields 
numerical QNM frequencies within less than 2% of the values in Table [llj Likewise, we obtain a numerical power-law 
tail of the form for the monopole, with p = —3.08 for a/M — and p = —3.07 for a/AI = 0.99 in good agreement 



with the prediction p = —3 obtained from the low- frequency expansion of the wave equation underlying Eq. (40). 

Additionally, we have performed a convergence analysis for the more challenging of these two cases, that of a 
highly rotating BH background with a/M = 0.99. This analysis provides the estimate of the numerical uncertainties 
A^'ii/^'ii < 8% for the I = m = 1 mode at late times of the evolution and A^qq/^'oo < 3% for ^ = to = 0. 
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FIG. 2. Amplitude of the scalar field extracted at rex = 10 M as function of time obtained from the evolution of a scalar field 
initialized as spherical shell with space dependent mass = — lOM^/r^, around a non-rotating BH. 




FIG. 3. Evolution of a Gaussian profile of a massless scalar field with width w = 2 M centered at ro = 12 M around a 
Schwarzschild BH (left panel) and around a Kerr BH with a/M — 0.99 (right panel). We depict the / = (solid black line) 
and I = m = 1 (red dashed line) multipoles. 



C. Massive scalar fields 



Having tested the numerical framework, we next explore the dynamics of generic, massive scalar fields. A mass 
term introduces a new scale to the problem and new features appear in the evolution which depend on the particular 
details of the initial configuration. Roughly, these configurations can be classified into three groups. 



Bound state configurations. These configurations are characterized by unusually long-lived modes [211 1331 - 



inZl SH IS21 IM] J and exist for any m, ^5 7^ 0. These modes are described by Eq. (38) for small M^s- If they. 



furthermore, satisfy the condition ujr. < i^s < mVLH they are subject to the supcrradiant instability. 



• Rapidly damped configurations. Generic initial profiles of massless or massive scalar fields typically decay on 
short time scales via quasi-normal ringdown followed by a power-law tail. 
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• Beating regime. Additionally, massive fields may exhibit long-lived, strongly modulated oscillations which result 
from the interplay between different overtones of the same multipole. 

In the following we will discuss numerical evolutions of initial configurations for each of these classes in more detail. 



Bound states 



- fundamental mods 

- overtone mode 





FIG. 4. Left and center: Initial profile of I'J'iip for a bound state m — 1 dipole configuration with Af/ig = 0.42 in a 
Schwarzschild (left panel) and a/M — 0.99 Kerr background (mid panel). Black solid lines represent the fundamental mode 
and red dashed lines the first overtone. Right: Relative change of the modulus of the (1, 1) mode extracted at r^x = 20 M. 



We have constructed bound-state initial configurations for fields with a mass coupling Af/15 = 0.42, cf. Eq. (30 1, 
following [9l [Ml |32], and evolved them in a Schwarzschild or a Kerr background with a/M = 0.99. In the left and 
center panels of Fig. [4] we show the modulus of the fundamental mode and first overtone of the m — 1 dipole bound 
state along the a;-axis; cf. Table |ll] For non-rotating BHs the fundamental bound state is localized near the origin, 
whereas its maxima are shifted to larger radii as the rotation parameter a/M is increased. For a highly spinning BH 
with a/M = 0.99, the fundamental mode is peaked at around r ^ 12 M. For the first overtone we observe a node 
at ?'nodo ~ 22.5 M and at rnode ~ 26.5 M for a/M = and a/M = 0.99, respectively. The specific structure of the 
bound-state profiles will play an important role in our analysis of beating effects below. 

Throughout the time evolution, we expect the bound state scalar field to remain localized in the vicinity of the BH. 
Thus, by construction, its absolute value ^'n^'^j^ ~ exp(— iwt) exp(zaji) const, should remain almost constant in 
time, with a small growth rate of Mojj ~ 1.5 • 10"'' |361 137j . This behavior is confirmed in the animations generated 
from our numerical data and made available online [SS]. Here we have tested these properties numerically by extracting 
the dipole mode as a function of time at rex = 20 M. The result is shown in the right panel of Fig. [4] for a Kerr 
background with a/M — 0.99. The scalar field varies by less than ^ 2% until t ^ 200 M and by less than ^ 1% at 
late times, which is within the numerical uncertainties. 

We note that bound states are unstable states, but have a long instability time scale of ~ lO^M (see Table [ll| , 
about three orders of magnitude larger than the evolution times feasible within our framework. Over the time range 
covered in the figure the instability has not yet generated a visible growth in amplitude. 



2. Damped states: ringdown and tails 



In order to study the behaviour of rapidly damped configurations, we initialize the field by a Gaussian wave pulse 
with ro = 12 M and w — 2 M according to Eq. 



(24a 



The specific choices of the mass coupling Miig, rotation rate 
a/M and the initial mode contributions for our set of simulations are summarized in Table |l] The resulting dipole 
and monopole amplitudes for a subset of our simulations are shown in Fig. [5] For all simulations we observe the 
expected pattern of an early transient followed by quasi-normal ringdown and a late-time tail which is dominated 
by an oscillatory behaviour characterized by the mass term according to Eq. ([39| and governed by scattering off 



spacetime curvature according to Eq. (41b I at late times 



At intermediate times we find for the case of an to = 1 dipole with Mfis = 0.1 (left upper panel of Fig. [5]), an 
oscillatory decay of the field as ~ ^-2-543 gin(O.lt), within 2% of the expected value for tails at intermediate times 
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FIG. 5. Upper left: I = m — 1 multipole of run sS_m001 extracted at lOM. The solid (black) line refers to numerical data, 
the dashed (red) to an oscillatory tail fit. Upper right: Same for sS_mlOO, extracted at Tcx ~ 20 M (black solid line) together 
with the tail fit (red dashed line, we omit the oscillatory term for clarity). Lower Panels: I — m — multipole of runs sS_m042 
(left, extracted at 20M) and sK_m042 (right, extracted at 25M). The solid black lines refer to the numerical data while the 
red dashed lines denote the the fit to the envelope of the oscillatory late-time tail. 



; cf. Eqs. (39 1 and (41a). At very late times, we expect a tail of the form (39) with p = —5/6 to dominate, but 



simulation of sufficient duration are computationally too expensive with our present computational framework. 



Furthermore, we consider larger mass couplings Mfis — 0.42 and Mfis — f. We present the I = m = 1 dipole 
mode of sS_mlOO in the top right panel of Fig. [5] and the monopole modes for both runs sS_m042 and sK_m042 in 
the bottom panels of Fig. [s] Here, the intermediately-late time tail with p — —{I + 3/2) in Eq. (39) appears to be 
suppressed. Instead, the decay with p — —5/6 expected at very late times clearly dominates the signal. Focusing 
on the case with mass coupling M/xg — 0.42 we observe that the tail is present for both spinning and non-spinning 
BH backgrounds; cf. bottom panel of Fig. [5] For a/M = and a/M = 0.99, we find the exponents p — —0.880 and 
p = —0.861, respectively, which agree within 5% with the theoretically expected late-time behaviour j5TH85] . We 
note that these results are independent of the extraction radii r^^ at which the field is observed. 



15 




1000 2000 








■0 2000 4000 6000 8000 10000 



2000 4000 6000 8000 10000 

t/M 



FIG. 6. Real part of the m = 1 dipole obtained at selected extraction radii rex for a massive scalar field with Mjis = 0.42 in 
a Schwarzschild (upper panels) and Kerr background with a = 0.99M (bottom panels). The location of the node of the first 
overtone is Tcx = 22.5 M for the Schwarzschild case (upper left most panel) and rex = 26.5 M for the Kerr case (bottom center 
panel), rex ~ 50 M corresponds to the overtone's local maximum in the Kerr case. 



3. Massive scalar fields: mode excitation and beating 



In addition to the well-known ringdown and decay, the evolution of massive Gaussian wave packets around BHs 
can exhibit more complex patterns. In particular, we expect massive scalar fields in Kerr backgrounds to eventually 
show exponential growth if they satisfy the superradiance condition Eq. ([!]) as well as wjj < ^g. The timescale for 
this instability, however, is of the order > 10^ M which is computationally too expensive to be realized within our 
numerical framework. In the following, we therefore focus on the complex signals observed at earlier times up to 
^10* M in the evolution of massive scalar fields around Schwarzschild or Kerr black holes. 

In particular, the time evolution of the I = m — 1 mode of the scalar field with mass coupling M/ig = 0.42 exhibits 
an oscillatory pattern with significant modulation of the amplitude. The quantitative behaviour of the oscillations, 
however, depends sensitively not only on the initial data but also on the radius where the modes are measured. This 
is illustrated in Fig. [6j where we show the m = 1 dipole extracted at different radii for Mfig = 0.42 considering a 
Schwarzschild or a Kerr BH background with a/M — 0.99. For each configuration shown in the figure, we have chosen 
three extraction radii, including the location Tnodc of the node, 22.5 M and 26.5 M respectively, for the Schwarzschild 
and Kerr case. 

This phenomenon can be explained in terms of a beating modulation between two or more long-lived modes 
described by Eq. (38 1. This beating modulation depends on the relative strength of the different overtones which, in 



turn, depends on the extraction radius; for the case of a guitar string, for example, a specific mode cannot be excited 
at the location of its nodes. Our choice of initial parameters for this example did not involve any finetuning and we 
expect these features to be present in the time evolution of generic massive fields around BHs provided only that at 
least two long-lived modes are excited. 

Our interpretation is confirmed by the Fourier spectra of the m — 1 dipole obtained at different radii which are 
shown in Fig. [t] for the two background spacetimes. For the Kerr case with a/A/ — 0.99 (right panel), we see that 
the signal is dominated by the fundamental mode at r^x — 20 M < rnodo and has two overtones of smaller amplitude. 
This is in agreement with the weakly modulated high frequency signal in the time evolution in the bottom left panel 
of Fig. [6] As expected, the amplitude modulation is particularly weak for rox = ''node = 26.5 M which coincides with 
the node of the first overtone. In fact, the small amount of modulation visible in the bottom center panel of Fig. [6] 
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FIG. 7. Spectra of the I = m = 1 mode of the massive scalar field with M^s = 0.42 evolved in the background of a 
Schwarzschild (left panel) or Kerr BH with a/M — 0.99 (right panel). The lines correspond to the waveforms measured at 
different extraction radii. In particular, Vex = 22.5 M (left) and r^^ = 26.5 M (right) correspond to the nodes of the first 
overtone. 



at early times is mostly due to the second overtone which, however, damps out at later times. Then, the envelope 
is almost constant with a small growth rate Mujj 10~^, in order-of-magnitude agreement with previous studies 
[37l|44]. The situation is markedly different at = 50 M, where the first overtone has a local maximum and the 
comparable strength of fundamental mode and overtone result in the strong beating modulation displayed in the 
bottom right panel of Fig. [6j This analysis can be repeated for the non-rotating case, with similar conclusions. In 
contrast to the Kerr case, however, all modes decay resulting in overall damped signals. 

Before proceeding with a detailed analysis of the beating modulation and mode excitation, we perform a convergence 
analysis of this setup which represents the most demanding and longest of our simulations. For this purpose, we have 
evolved the setup sK_m042 using three different resolutions he — M/QO, h„i — M/72 and hh = Af/84. We present 
the convergence plot, i.e. the differences between the coarse-medium and medium-high resolution runs in Fig. [8] for 
the monopole (left panel) and dipole (right panel) . Our numerical results show second order convergence throughout 
the simulation. We estimate the discretization error to be about Aipn/tfin < 1.1% for the entire interval and 
^V'oo/V'oo < 1% I = m ~ mode at early time which increases to AipQo/^oo ^ 7% at late times. 
Beating. In order to better understand the beating pattern quantitatively, let us for simplicity consider the presence 
of only two modes with similar frequencies: the long-lived fundamental mode with frequency uiq = utR^ + iujj q ~ (jJRfi 
and amplitude Aq and the first overtone with frequency ui = ujr,o + 5io, where Sio <^ 1, and amplitude Ai. The 
generalization to a larger number of modes is straightforward. For illustration, we explicitly list dio for a selected 
subset of configurations in Table |llj Because all these modes are long-lived, they are well approximated by pure 
sinusoids over short periods of time. The waveform ^' is then described by the superposition 

^ r^{Aa- Ai) sin(wfl,oi) + sin{ujR^ot) cos(5iot/2) . (42) 

The modulation in amplitude is governed by the low-frequency signal cos SiqI whereas the beating amplitude depends 
on the relative strength of the modes Aq, Ai. For equal amplitudes Aq — Ai, for instance, the total signal is given 
by a sinusoid modulated by a cosine. The beating frequency Siq can be estimated by fitting the envelope of the total 
signal. For the two cases displayed in Fig. [6] we obtain Siq ^ 0.0074 (for sS_m042) and Sno ~ 0.0063 (for sK_m042), 
in excellent agreement with the corresponding predictions listed in the third and fourth row of Table |ll] 

This picture is confirmed by calculating the mode spectra from a Fourier transformation of the time series. The 
spectra thus obtained for the two configurations and different extraction radii are shown in Fig. [7] and each exhibit 
three pronounced peaks which correspond to the real parts of the fundamental and first two overtone QN frequencies. 
A spectral analysis applied to the Kerr case a/M = 0.99 reveals Siq = 0.0063 and S20 — 0.0087. Thus, the evolution of 
Gaussian initial data excites a third long-lived mode which is not listed in Table |TTj We have verified that this mode 
indeed exists, via a continued fraction scheme in the frequency domain, corresponding to Muj = 0.41730-|-i2.265x 10^^. 
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FIG. 8. Lower panels: Convergence plot of the I = m — (left panel) and I — m = 1 (right panel) modes of simulation 
sK_042. We present the differences between the coarse-medium and medium-high resolution runs, where the latter is amplified 
by Qi ~ 1.66 indicating second order convergence. 



Mode excitation. The mode excitation can be put on a more rigorous framework: Leaver's seminal work, in 
particular, has estabUshed important results in this context and for further details we refer the reader to the original 
work [13] as well as comprehensive follow-up studies [HI [Ml [HZ]- The upshot is that each quasi-normal mode, which 
corresponds to a pole in the complex-frequency plane, is excited to a different degree depending on the initial data 
and on the mode in question. The QNM contribution can be isolated from other features of the signal, such as 
the late-time tail, using the Green's function technique [MllSHllH^. In this formalism, the scalar field amplitude at 
intermediate times is given by a sum over quasinormal modes as 

* = ^C„e-'""*VnK,r), (43) 

where '0„(w„, r) is the quasinormal mode eigcnfunction and aj„ its frequency, both quantities can be computed via the 
Fourier-domain ordinary differential equation that governs massive scalars in the Kerr background [141 1861 187] . The 
numbers C„, called excitation coefficients characterize the amplitude to which each mode is excited. Two quantities 
are crucial to determine the excitation coefficients [87]: the behavior of V'n close to the eigenfrequency w„, and the 
convolution of the eigenfunction with the initial data. Thus, for instance, the relative amplitude between different 
modes depends strongly on the point where this amplitude is evaluated: if it is close to a node of one of the modes, 
the mode in question will have a very small amplitude: by definition a mode is not excited at its node. Likewise, 
localized initial data close to the node of the mode do not excite the mode in question, a well-known result for closed 
systems [57] . 

We have not attempted a complete quantitative understanding of mode excitation for this work, a preliminary 
analysis indicates that the excitation coefficients are indeed of comparable magnitude away from the nodes of the 
eigenmodes, but vary substantially close to the nodes. 



IV. PROCA FIELD EVOLUTIONS 



(8a 



The evolution of vector fields, as for example the photon, are governed by Eq. (8b I which, after substitution of the 
V^Ay — \7,yAfj_ and the Lorenz condition, Eq. (14), becomes rather similar to its scalar counterpart 
For this reason, it has been believed for a long time that massive vector fields should also be prone to a "BH 



definition F,, 



bomb"-like superradiant instability. In fact, Rosa and Dolan [32] conjectured that instabilities of vector fields should 
not only be present but can be much stronger than those of scalar fields. This has recently been verified explicitly by 
Pani et al [40l [41] who derived the instability growth rates of massive vector fields in the slow-rotation approximation 
and found them to be several orders of magnitude larger than their scalar field counterparts. Their results, however, 
have been obtained only in the regime of slowly rotating black holes. Even though their results are conjectured to 
hold for arbitrary spins, a definitive answer to this question calls for the modelling of vector fields in generic Kerr 
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Run 


a/M Mfiv 




) w/M 


Grid Setup 


vlS_mOOO 


0.00 


0.00 


Supl 


2.0 


{(192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 


v2S_m000 


0.00 


0.00 


Supl 


30.0 


{(192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 


vlS_m010 


0.00 


0.10 


Supl 


2.0 


{(192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 


vlS_m020 


0.00 


0.20 


Supl 


2.0 


{(192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 


v2Kl_m040 


0.50 


0.40 


Sup2 


30.0 


{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 


vlK2_m000 


0.99 


0.00 


Supl 


2.0 


{(192, 96, 48, 24, 12, 6, 3, 1.5), h = M/96} 


v2K2_m040 


0.99 


0.40 


Supl 


30.0 


{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/64} 


v2K2_m042 


0.99 


0.42 


Supl 


30.0 


{ (1536, 384, 192, 96, 48, 24, 12, 6, 3,1.5), h = M/64} 


v2K2_m044 


0.99 


0.44 


Supl 


30.0 


{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/64} 


vlK2_ml00 


0.99 


1.00 


Supl 


2.0 


{(1536, 384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/60} 



TABLE III. Initial setup for simulations of Proca fields with mass coupling M^v in BH background with dimensionless spin 
parameter a/M . The initial Gaussian pulse with width w/M is located at ro = 12 M and consists of a superposition of s = — 1 
spin- weighted spherical harmonics given by Eq. \44a.\ ("Supl") or Eq. \44h\ ("Sup2"). We further denote the grid setup, in 
units of the BH mass M following the notation of Sec. II E in [72] . 



geometries. In the frequency domain, such modelling represents a formidable challenge because the equations of 
motion appear to be non-separable. Here we therefore address this question in the framework of numerical evolutions 
in the time domain, where the non-separability of the equations does not represent a serious obstacle. 

For this purpose, we prescribe initial data in the form of a superposition of Gaussian pulses centered around 
tq = 12 M and composed of several multipoles with angular dependence given by s = — 1 spin- weighted spherical 
harmonics -lYim- Specifically, we have chosen linear combinations 



( -lYi-i + -lYii) - ( -1F2-1 - -1Y21) - -lYio - 



1 
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1 



( -1Y2. 



15^22) 



V2 



( -1Y3-2 + -1132) , 



3n/5 ' 3V35 

( -1^1-1 + -i^ii) — ( -1Y2-1 — -1Y21) , 



(44a) 
(44b) 



such that the imaginary parts of the fields vanish. By virtue of the evolution equations (20a)- (20c), purely real initial 



data remain real throughout the evolution so that our choice reduces the computational requirements to evolving 
only seven - instead of 14 - independent variables. We have summarized the initial configurations of our set of 



simulations together with the grid setup in Table III In order to analyse the time evolutions of the vector field 



we decompose its time component </? and the Newman- Penrose scalar $21 constructed from the spatial components 
according to Eq. (32), into multipoles by projecting them onto spherical harmonics with spin weight s = and s = — 1, 
respectively. 



A. Massless vector fields 



We first study the behaviour of massless vector field perturbations with M/iy = 0. This case has been studied 
extensively in the literature [2| and therefore also enables us to compare our findings with previous investigations. 
Our results for massless Proca fields are summarized in Table IV and Figs. [O] and [Tol Let us first consider the simplest 
case of a massless vector field in Schwarzschild background. The time evolutions of the I ^ m = 1 multipole obtained 
for initial Gaussian pulses of width w = 2 M and w = 30 M are shown in Fig. [9j For the narrow pulse (solid curve 



in the figure) we clearly identify the pattern familiar from our scalar field evolutions in Sec. IIIB an early transient 
whose details depend on the initial configuration is followed by a quasi-normal ringing characterized entirely by the 
BH parameters and a late-time tail. For an initially broad pulse (red dashed curve), however, the initial transient is 
directly followed by a power-law tail with no visible intermediate ringdown stage. This feature has been reported for 
scalar fields in Refs. [SSI US] and is a consequence of the negligible excitation of the long-lived fundamental mode and 
low overtones by broad pulses; the high overtones which are excited significantly by this type of initial data rapidly 
decay before the transient gives way for a clear QN ringdown pattern to emerge. 

For a quantitative comparison with calculations performed in the frequency domain, we have fitted the ringdown 
part of our numerically extracted multipoles for the case of a narrow initial Gaussian with exponentially damped 
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a/M ilm) 



0.00 (10) 0.2483 - 10.0925 0.248 - j0.092 
0.00 (20) 0.4576 - j0.0950 0.455 - j0.093 



0.99 


(10) 0.2743 - 


10.0759 


0.274 - 


zO.075 


0.99 


(11) 0.4634- 


jO.0313 


0.464 - 


20.033 


0.99 


(20) 0.4999 - 


zO.0800 


0.498 - 


i0.079 


0.99 


(22) 0.9099 - 


10.0301 


0.887- 


10.037 



TABLE IV. Quasi-normal mode frequencies of massless vector perturbations in a Schwarzschild or Kerr BH background with 
a/M = 0.99. The values for Mui^^ have been obtained from fits to our numerical evolution of the field, whereas those for 
Mojif^ have been computed with the continued fraction method [9l 1111 188] . 



sinusoids. The resulting complex frequencies W;™™ are listed in Tabic 
from frequency-domain calculations [9l [TTl [Ml [88] . 



IV 



and agree well with the values obtained 
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FIG. 9. The I — l,m — multipole of $2, extracted at rex ~ 10 AI from the time evolution of a Gaussian pulse of width 
w = 2 M (black solid line) and w — 30 M (red dashed line) around a Schwarzschild BH. 



The time evolution of a massless vector field initialized as a narrow Gaussian of width w = 2 M around a rapidly 
spinning Kerr BH with a/M = 0.99 is displayed in Fig. 10 and qualitatively agrees with the corresponding simulation 
around a Schwarzschild background. Note, however that the quasi-normal ringdown is significantly slower for the 
spinning case which is also reflected by the relatively small imaginary quasi-normal mode frequencies listed for this 
case in Table |IV[ Accurately measuring such small imaginary components in numerical simulations represents a 
considerable challenge which is why we have chosen a higher numerical resolution for this particular simulation; 
cf. Table III We thus obtain agreement of a few % with frequency-domain predictions for the imaginary part. In 
contrast the real part of the frequency is much easier to extract numerically and shows excellent agreement around 
1 % or less with the values calculated in the frequency domain. 



B. Proca field in Schwarzschild backgrounds and in the slow- rotating limit 



We now consider time evolutions of massive vector, or Proca, fields in a Schwarzschild background. For this purpose, 
we choose the mass parameters Mfiy = 0.1 and M/iv = 0.2 which allow for a direct comparison with recent QNM 
computations by Rosa & Dolan [42 and Pani et al [401 [41]. In contrast to the massive scalar field, the Proca field has 
three degrees of freedom and the resulting radiation multipoles can be classified into three groups: axial modes with 
spin s = and two polar modes with s — ±1. Following Refs. [40 42], we denote s = multipoles as even scalars (ES) 
and the s = -1-1, —1 as odd (O) and even vector (EV), respectively We emphasize that these three different degrees 
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FIG. 10. Time evolution of the I — m = 1 mode of $2 of run vlK2_m000, i.e., a massless vector field in a Kerr background 
with a/M = 0.99 extracted at rox = 10 M. 





(ES) Mujfo 


(0) M4g 


(EV) Mu^Yr i^i) 






P {v>io) P (*2,;o) 


0.1216 - 


10.0791 




0.125- 


10.080 






1.51 


1 0.3054 - 


10.0914 0.2435 - 


10.0944 0.2539 - 


10.0887 0.245 - 


zO.098 


0.252 - 


i0.087 


0.49 0.49 


2 0.4921 - 


10.0943 0.4552 - 


jO.0955 0.4610 - 


jO.0938 0.462 - 


10.091 


0.452 - 


i0.096 


0.48 0.49 



TABLE V. QNM frequencies and tail exponents for Proca field perturbations in the Schwarzschild background, for Mfj.v = 0.1. 
Modes b/'^ have been computed with the continued fraction and forward-integration method in the frequency domain |40H42j . 
These modes are divided in odd (O), even scalar (ES) and even vector (EV) parity; see text for details. The QNM frequencies 
extracted from our numerical simulations are shown as ct)""™, as extracted from ipiQ and $2,;o. Finally, we also show the decay 



exponent p of the oscillatory tail; c. f. Eq. ( 45 1 



of freedom have different spectra except for the massless case where the spectra of the vector modes are degenerate 
while the scalar mode reduces to a gauge field [21 SH HSj . The quasi-normal ringdown frequencies oof^ obtained from 
frequency-domain calculations for the three types of multipoles are shown in Table [V| for M/xy = 0.1 considering a 
Schwarzschild background. 

In order to obtain numerical estimates for the frequencies, we have evolved Gaussian initial data of width w — 2 M 
centered around ro = 12 M for our two choices of M/iy =0.1 and 0.2. The resulting dipoles of $2 extracted at 
Tex = 10 M are shown in Fig. [TT] and reveal the by now familiar pattern of early transient, ringdown and tail. Fitting 
an exponentially damped sinusoid to the ringdown part of the dipoles of $2 (solid curve in the figure) as well as the 
scalar component ip for the case M/i = 0.1 yields numerical estimates for the frequencies listed in Table [V] as W;™™. 
These estimates agree with the frequency-domain predictions within a few percent or less. Note, however, that the 
frequencies for some types of modes are very similar, so that we cannot unambiguously identify which modes are 
excited by our particular choice of initial data. For I = 2, for instance, our numerical results are compatible with 
both, an odd or even vector mode. 



In Fig. 11 we see that from time t ^ 100 AI onwards, the signal is dominated by the tails whose functional form is 



given by a decaying sinusoid of the form [90] 

$ ^ t-C+p) sin(^yt) , (45) 

at intermediate times, where p = 3/2 + s depends on the spin, or parity, of the mode. At late times, on the other 
hand, the signal is expected to follow the universal behaviour 

$ - sinipvt) , (46) 

independent of the spin s of the field '5P, [90] . Numerical estimates for the exponent p extracted from our numerical 
results for $2 and (p are shown in Table [v| and are in excellent agreement with the prediction p = 3/2+ s for intermediate 
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times and spin values s = for the monopole and s = —1 for dipole and quadrupole. A similar analysis for the larger 
mass parameter M/z = 0.2 also leads to good agreement between numerical and frequency-domain results, albeit with 
slightly larger discrepancies of < 7 %. 




FIG. 11. Time evolution of the Proca field with mass coupling Mjiy = 0.1 (black solid line) and Mfiy = 0.2 (red dashed 
line) in a Schwarzschild background. We show the I = l,m = multipole of the Newman- Penrose scalar $2, extracted at 
r,, = 10 M. 



Before we discuss in detail the behaviour of Proca fields in rapidly rotating Kerr backgrounds, we briefly test our 
code in the case of a smaller rotation rate a/M — 0.5 where the slow- rotation approximation of [40l|41] is expected to 
provide rather accurate results. For this purpose, we have evolved a Proca field with a mass parameter M^y = 0-4 
(simulation v2Kl_m040 in Table III I which can be shown to result in a stable mode, i.e., it does not satisfy the 



superradiance criterion ([T]). For this configuration we observe an extended transient resulting from the wave packet 
impinging on the black hole followed by a slowly decaying QN ringdown signal after t ~ 200 M. Fitting a damped 
sinusoid to this ringdown part, we find a ringdown frequency Mlou = 0.389 — i0.0023 in excellent agreement with 
semi-analytic calculations in the slow-rotation approximation [401 141| . 



C. Instability of Proca fields in highly spinning Kerr backgrounds 

According to the superradiant condition ([T]), rapidly rotating black holes are most likely to induce superradiance 
phenomena in ambient vector fields and we therefore discuss the case of Proca fields in a Kerr background with 
a/M = 0.99 in this section. In fact, recent calculations [40, .41^ indicate that the maximum instability growth rate is 
realized in this background spacetime for a mass parameter of about Mfiy = 0-4 and our discussion will focus on this 
case supplemented by additional simulations with M/zy = 0.42, 0.44 and 1.0. The time evolution of the I — m — 1 
multipole of the Newman-Penrose scalar $2 as well as the scalar component ip obtained for the case Mfiy — 0.4 is 
shown for several choices of the extraction radius in Fig. |12| Animations of the numerical evolutions are available 
online [85 . 

Beating of modes. The time evolutions of the Proca field shown in Fig. [12] reveal strong amplitude modulations 
similar to those observed for scalar fields in Fig. [6j They also depend sensitively on the location of the measurement 
specified by the extraction radius rex and, again, we interprete this feature as a beating effect. In order to study 
this in more detail, we have computed the Fourier spectra of the waveforms at different radii and show the results 



in Fig. 13 The curves clearly reveal several peaks corresponding to separate mode contributions (fundamental mode 
or overtones) of the dipole field and the relative amplitude of these peaks varies significantly with extraction radius. 
Consider, for example, the spectrum of the dipole of $2 in the left panel of Fig. [TSj At — 10 M the amplitude 
of the lowest frequency mode is about two thirds of the amplitude of the strongest peak near Mujr = 0.39, but its 
relative strength rapidly decreases at larger extraction radii. This is also refiected in the observed time evolutions 
of the $2 dipole in the upper panels of Fig. [12] the waveform extracted at Tex = 10 M exhibits a high frequency 
modulation which is weakly present at Tqx = 25 M and entirely absent in the waveform measured at Tqx = 40 M . 
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Likewise, the high-frequency modulation of the time component Lp weakens at larger radii as the relative amplitude 
of the two lowest peaks in the spectrum (right panel of Fig. 13 1 decreases. 

There remains one important issue: which conditions lead to the beating of modes? While beating is expected to be a 
generic feature, it is not always excited. Here, we attempt to determine indicators for the (non-)observability of beating 
effects. Because we have performed a multipole expansion of the fields, and beating is present on a fixed multipole, 
then clearly a beating pattern requires the excitation of different overtones. For this reason, a modulation of the 
amplitude is only triggered by evolutions of generic initial pulses |^ Furthermore, our numerical simulations indicate 
no dependence on the background spacetime; Fig. [6] shows beating effects of the scalar field around Schwarzschild and 
Kerr BHs. Beating is therefore not related to superradiance. On the other hand, the excitation of different modes 
is merely a necessary condition but not sufficient to trigger amplitude modulations. This becomes evident in the 
evolutions of Gaussian initial data for the vector field with small mass parameters in a Schwarzschild background 
which only exhibit a quasi-normal ringdown and tail. In order to shed further light on this question, we therefore 
consider the entire set of Proca evolutions in BH backgrounds and summarize our results as follows: (i) We do not 
observe beating nor any long-lived modes for mass parameters A//iy = 0.1,0.2 in Schwarzschild background; (ii) we 
observe beating in the dipole mode for M/iv = 0.4, 0.42, 0.44 in rapidly rotating Kerr; (iii) we observe beating in the 
quadrupole but not in the dipole for M/iy = 1.0 in the background of a highly spinning Kerr BH. Let us first consider 
in our interpretation of these observations the case of small mass parameters M^sy ■ The absence of long-lived modes 
in our simulations is most likely a consequence of low-mass bound states being concentrated far from the black holes. 
For our choice of initial Gaussian pulses centered around rg = 12 M, these states are therefore only weakly excited and 
play no significant role in the evolution. Furthermore, beating patterns for small mass parameters have particularly 
large periods. We can therefore not rule out that amplitude modulations might become visible in these evolutions 
at late times far larger than the evolution times feasible with our numerical framework. We suspect a similar reason 
is behind our not observing beating modulation in the dipole mode in case of large mass coupling M/iy = 1. We 
therefore tentatively conclude that two conditions need to be met in addition to the presence of at least two modes 
in order to observe beating on time scales of ^ 10"* M . (i) the bound states must not be localized far away from the 
peak of the initial data and (ii) beating periods need to be sufficiently short. 

Instability of Proca fields around Kerr. The most striking feature of the time evolutions in Fig. [12] is the 
amplitude growth of the signal over time (ignoring the early transient stage) . We interprete this growth as a signature 
of the superradiant instability of massive vector fields and estimate the growth rate of the m = 1 dipole of $2 obtained 
for Miiv = 0.4 and a/M = 0.99 as 

Mw/ - (5 ± 1) • 10"'' ^ ^ - (3.3± 1.7) X 10^ (47) 

and about half that value for the instability rate of (/?. We believe the discrepancy of the growth rates of $2 and ip 
is due to the different growth scales of the axial and polar sector and the fact that Lp only sees the axial sector. We 
note that our estimates for the growth rates are of the same order of magnitude as the value r/M ~ 10'^ derived 
from extrapolation of slow-rotation calculations; cf. Eq. (98) in Ref. [H] and Fig. 2 in Ref. [ID]. The superradiant 
instability time scales for vector fields are thus up to four orders of magnitude larger than those of their scalar counter 
parts which renders possible their identification in the numerical evolutions presented in this work. 



V. CONCLUSIONS 



We have performed an extended study of the time evolution of massless and massive scalar and vector fields in 
Schwarzschild and Kerr BH background spacetimes with spin parameters up to a/M = 0.99. Our results are consistent 
with published results obtained in the frequency domain in so far as these are available. For evolutions involving 
exclusively short-lived modes, we observe the known pattern of an initial transient followed by an exponentially 
damped ringdown phase and late-time tails. 

For the case of massive fields, there exists a second class of long-lived modes called bound states whose evolution 
exhibits a significantly richer structure. In particular, we have identified a beating pattern caused by the interference 
of different modes (the fundamental mode and overtones) as manifest in the Fourier spectra calculated from the time 
evolutions. The relative amplitude of different modes in the spectra and, thus, the specific shape of the beating 
pattern strongly depends on the observation radius. We believe that these beating effects provide an explanation 
for an apparent discrepancy between time and frequency domain calculations of instability growth rates of scalar 



as opposed to pure bound states. It is important to note that by pure bound states we mean stationary solutions of the linearized field 
equations in the Kerr background; therefore effects such as mode coupling are already taken into account: pure modes do not couple to 
other multipoles when the expansion basis is taken to be spheroidal harmonics. 
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FIG. 12. Time evolution of the Proca field with mass couphng Mfj,v = 0.40 in Kerr background with a/M = 0.99, at different 
extraction radii. We plot the I = m = 1 mode of the Newman- Penrose scalar $2 (upper panels) and of the scalar component 
ip (lower panels). 
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FIG. 13. Spectra of the I — m = 1 mode of the Newman-Penrose scalar $2 (left panel) and of the scalar component ip 
(right panel) excited by the Proca field with mass coupling Mfiv = 0.40 in Kerr background with a/AI = 0.99. Different lines 
correspond to the waveforms measured at different extraction radii. 



fields around Kerr backgrounds. Specifically, Ref. [60] (see their Fig. 3) reported a time evolution lasting about 
3 X lO'^M of a massive scalar field with Mf^s = 0-25 in a Kerr background with a/M — 0.9999 and obtained a growth 
timescale approximately two orders of magnitude smaller than that predicted by frequency-domain calculations. For 
this particular configuration, the frequency-domain estimate, Eq. (38 1, predicts a beating period r ^ 5 x lO'^M, about 
the entire duration of their time evolution, and we believe that their estimate of the growth rate was correspondingly 
contaminated by the resulting amplitude modulation. 

Based on our numerical results, we conjecture that two conditions are required for the presence of beating: (i) 
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generic initial configurations of massive fields as opposed to single-frequency bound state initial data, and (ii) initial 
data able to excite these bound states. For extremely small masses for instance, the bound states are localized far 
from the BH. Accordingly, initial data peaked far from the BH would be required to excite them. 

In an accompanying paper Dolan [61' has investigated the massive scalar field instability as well as a system 
involving massless scalars but enclosed by a mirror in the time domain. Specifically, the author employs a "coupled 
1 + 1 dimensional" numerical scheme which allows for significantly longer evolutions times up to t ~ lO^M. This work 
finds bound state frequencies and instability time scales in excellent agreement with frequency domain calculations 
and, due to the generic initial configurations, confirms/ supports our findings of the beating phenomena. 

Finally, we have confirmed the existence of superradiant unstable massive vector fields around rapidly spinning 
BHs. In contrast to massive scalar fields, the instability of Proca fields is stronger by about four orders of magnitude 
leading to growth times as short as t/M ~ 3.3 x 10^ in rapidly spinning Kerr BH backgrounds. It is instructive 
to translate this number for the case of realistic astorphysical black hole candidates. For a solar mass BH and a 
supermassive BH of the size of SagittariusA* [M - 4.1 • lO^M©) at the center of our galaxy we obtain timescales of 
T ^ 9 ms and t ^ A - 10'* s, respectively. 

It has been argued by Pani et al. [40 that the angular momentum thus extracted from the black hole provides a 
mechanism to observationally constrain the mass of the photon. In view of the short time scales even for a supermassive 
black holes, spin measurements of the BH at the center of the Milky Way, may indeed constrain the photon mass to 
unprecedented levels. The exciting prospect of using large, supermassive black holes to understand the microscopic 
world raises several questions: 

(1) Influence of accretion disks. Astrophysical black holes are not isolated, but typically are surrounded by accretion 
disks. Can the interaction with matter kill the instability? It was argued previously [40j that the superradiant 
instability is a global mode, on scales larger than the BH. On these scales matter is electrically neutral, and the 
coupling should be negligible. Therefore, it is not likely that (thin) accretion disks can quench the instability; thin 
disks are expected to lie along the equator of a Kerr BH and are not expected to interact strongly with boson clouds 
that extend well off the equator. The influence of thick disks or other effects is unknown at the moment, but it is 
surely important to study these effects in more detail. 

(2) Self-interacting scalar flelds: One class of interesting problems involves massive scalar fields whose dynamics are 
described by additional non-linear terms, modelling their self-interaction. This open issue has first been addressed by 
Yoshino & Kodama [H] who modelled the collapse of a so-called bosenova. 

(3) Backreaction effects: As far as we are aware, all studies exploring massive flelds have been performed in the linear 
regime. Therefore it is of utmost interest to explore the fully non-linear regime, which allows for the investigation 
of the backreaction of the spacetinie, such as the spin-down of the BH due to (subsequent) superradiant scattering. 
This type of studies would enable us to glance at the end-state of the superradiant instability or, possibly, equilibrium 
conflgurations. 

The present study marks crucial first steps to explore an entire playground of exciting future applications of massive 
fields in BH spacetimes. 

Animations of the evolutions can be found online [85] . 



VI. ACKNOWLEDGEMENTS 

We thank Sam Dolan, Paolo Pani and Joao Rosa for useful comments and discussions and for making their data 
available to us. We are indebted to Hideo Kodama and Hirotaka Yoshino for explanations regarding time-evolution 
of bound states. We thank Gaurav Khanna for useful correspondence. We are especially indebted to Sergio Almeida 
for all his hard work on the "Baltasar Sete-Sois" cluster. We thank all participants of the YITP-T- 11-08 work- 
shop on "Recent advances in numerical and analytical methods for black hole dynamics" for useful discussions. 
We thank the Yukawa Institute for Theoretical Physics at Kyoto University for their kind hospitality during the 
early stages of this work. This work was supported by the DyBHo~256667 ERG Starting Grant, the NRHEP-295189 
FP7-PEOPLE-2011-IRSES Grant, the C5ffS0-g95^i^ FP7-PEOPLE-2011-GIG Grant, the ERC-2011-StG 279363- 
HiDGR ERG Starting Grant, and by FGT - Portugal through PTDC projects FIS/098025/2008, FIS/098032/2008, 
CTE-ST/098034/2008, GERN/FP/123593/2011. H.W. is funded by FGT through grant SFRH/BD/46061/2008. U.S. 
acknowledges support by the from the Ramon y Gajal Programme and Grant FIS2011-30145-C03-03 of the Ministry of 
Education and Science of Spain. A.I. was supported by JSPS Grant-in- Aid for Scientific Research Fund (G) 22540299 
and (A) 22244030. 

Gomputations were performed on the "Baltasar Sete-Sois" cluster at 1ST, the cane cluster in Poland through PRAGE 
DEGI-7 "Black hole dynamics in metric theories of gravity" , on MareNostrum in Barcelona through BSG grant AEGT- 
2012-2-0005, on Altamira in Gantabria through BSG grant AEGT-2012-3-0012, on Gaesaraugusta in Zaragoza through 
BSG grants AEGT-2012-2-0014 and AEGT-2012-3-0011, XSEDE clusters SDSG Trestles and NIGS Kraken through 



25 



NSF Grant No. PHY-090003, Finis Terrae through Grant GESGA-IGTS-234 and the COSMOS supercomputer, part 
of the DiRAC HPC Facihty which is funded by STFC and BIS. The authors thankfuUy acknowledge the computer 
resources, technical expertise and assistance provided by the Barcelona Supercomputing Centre — Centro Nacional de 
Supercomputacion and by Audrey Kaliazlin for computational support and technical advice with COSMOS. 



Appendix A: Flux formula 



From the Lagrangian £, 



we obtain, under the Lorenz condition (14), the equations of motion 

(V^V^-/i|)^-^^*i^''"^i^M.-V^'(*) = (A2) 

(V^V^ ~ My) - i?/^^ - 2fcaxion *F''^V^^ = -J^^> (A3) 

and the stress-energy tensor 

T^.=r;t + T*, (A4) 

T^u ■■= Ff^cF," - Ig^.F^pF^f' + i?yA,,Au - ^ii^A'^A^g^, , (A5) 

\ (V,,4'*V,* + V^.W,**) - ^.g^, (V^'J'*V'^* + 4*** + 2V{^)) . (A6) 

Note that the Chern-Simons term does not alter the stress-energy tensor. 

We are concerned with stationary axisymmetric spacetimes. Let us denote by t^^ and respectively, the stationary 
and the axisymmetric Killing field. The rigidity theorem states that there exists a Killing vector field which is 
normal to, hence null on, the horizon "H^ . Such a Killing vector field is given by the combination of and (j)^, 

X'' = t^ + , (A7) 

where ^Ih denotes the angular velocity of the Killing horizon of the black hole, and for the Kerr metric 



with r_|_ being the horizon radius = M + V — in the Boyer-Lindquist coordinates. 
Let us introduce the energy current by 

and consider two time-slices Si and E2 C {J^(Si) \ Si}, which intersect the event horizon "H^ at i = ii and t — t2, 
respectively. (One may view t as the Killing parameter of f^). Then the energy flux F that flows into the black hole 
in the time-interval I = {"H^ : ti < t < is given by 

F = I dS^ - / dS^ - / dA^n" 

= - / dA^x^J^ = / dS{x^eT^,) . (AlO) 

Here n'^ denotes the normal to 7^+, dA^ the volume element of the horizon H'^, and dS the area element of the 
cross-section H = D S, and (•) expresses the time average along the horizon. Superradiant scattering occurs when 
the energy flux going into the horizon becomes negative: 

F<0. (All) 
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We consider below the flux with respect to the scalar field and the vector field separately. For the scalar field, the 
mode decomposition is defined by 



where uj denotes the frequency and m the angular quantum number. It implies, in particular. 



(A12) 



(A13) 



Then, substituting (A6) into (AlO) and using the fact that Qfivf^X^ = on the horizon, we can immediately find the 
flux formula for ^P: 



JH 



and read off the superradiance condition 



< a; < mCln 



(A14) 



(A15) 



For the vector field, substituting (A5| into (AlO I, we have 



I dS {x^F^o.fF,") + III I d {x^'A^t'^A,) 

JH JH 



where we have again used gfiuX^f^ = oi^ the horizon. 
For the vector field, the mode decomposition is defined by 



Now, noting 



{£tA)^, = -iioAf^ , {£xA)f^ = -i{uj - mQ.)A 



(A16) 



and using Eq. (A17I and the Lorenz condition (14 1 we find 

X^F^af'F^" = g"^£xAa£tAfs + x''Af,V°'\/a,{t''A^) + divergence terms . 
On the horizon integral, we ignore the divergence terms and then obtain 

F^= [ dS {Reix^ F^^t- F,"* ))+ t^v [ dS {Reix^A^f^A:)) 

JH JH 

^io{uj~mn) f dS {\A\^) + fi^y f dS{\x''A^eA,\)- f d5 (Rc (x^A^V"V„(t''A:))) 

JH JH JH 



(A17) 
(A18) 
(A19) 



If we impose 



X^A^ = on •H^ 



we get the simple formula, similar to the scalar field case (A14) 



F^ ^ uj{uj ~ mil) / dS{\A\^) 

JH 



(A20) 



(A21) 



(A22) 



and the superradiance condition for the vector field is again given by (A15). 
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Appendix B: Spin-weighted spherical harmonics 



Here, we list the spin- weighted spherical harmonics up to / = 2 in spherical coordinates {0,4>} and Cartesian 
coordinates {x, y, z} given by Eq. (25 1. In case of spin- weight s = we obtain: 



1 = 0: 



J nf 



^00 



(Bl) 



/ = 1: 



^10 



Ali- 



ces 6 = 



3 z 
All r 



sm W cos ( 



-'1-1 



Stt 

^11 ; ^1-1 ~ ^11 



^10 = 0, 



3 X 
Stt r 



Y' - 
^11 — 



sm U sm ( 



Stt r ' 



(B2a) 

(B2b) 
(B2c) 



I = 2: 



^21 



Y'^ 

1 n' 



167r 



(3cos2 6i- 1) 



5 



■ cos f sm t/ cos ( 



2-1 



15 
32^ 

yi?. 



sin^cos(20) 



2 

167r \^ 

5 xz 
Stt r2 



y 



2-1 



-'21 



327r 

yR 



-"22 ' 



y^ 



y-* 



Y^ 

^2-2 



0, 



cos W sm f sm ( 



15 
32^ 



'sin(20) 



Stt r2 
' 15 
8^72 



-y^ 
^22 



(B3a) 

(B3b) 

(B3c) 
(B3d) 



where cos(20) = cos^ (f) ~ sin^ 0, sm{2(f>) = 2 cos sin (/), cos(3(?!)) = 4 cos'^ — 3 cos and sin(30) = 4 cos^ (/) sin — sin 0. 
We further summarize the s = — 1 spin- weighted spherical harmonics up to / = 2, where we have also defined 



I = 0: 



1^00—0) -1^00 — 



(B4) 



I = 1: 



vR 
-IJ^IO 



-in^ = 



167r 



cos cost 



1) 



3 
16^ 

8^ 



sin (/)(cos 9 — 1 
sin 6 = 



Stt r 



IGtt 



167r 



cos -I- COS ( 



3 x{z — r) 
IBtt rp 



3 y{z-r) 
IGtt rp 

_iy/o = 



3 x{z + r) 
IGtt rp 



sin 0(1 -I- COS 6*) = — 



3 y{z + r) 
IGtt rp 



(B5a) 
(B5b) 
(B5c) 
(B5d) 
(B5e) 
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I = 2: 



-1^2-2 — 
-1^2-2 = 
-1-^2-1 — 
-1^2-1 = 



-1^20 



-1J^21 
-1^21 



-1^22 



-1^22 - 



sin 6l(cos6'- l)(2cos2(^- 1) 



cos(/)sin(^sin6'(cos6' — 1) = — ., 
An V 47r r^p 



5 (z ~ r)(a;^ — y^) 
167r r^p 
5 a;?/(z — r) 



cos (/i(2 cos^ 6* — cos — 1) = 



5 x(2z'^ — zr — r^) 



IGtt 



j.2p 



IGtt 
15" 



■ sin <7!i(2 cos 6 — cos 6* — 1) = — 



5 y{2z^ — zr — r^) 



j.2p 



cos a sni ( 



15 zp 



167r 



IGtt 



cos </)(2 cos^ ^ + cos 6' - 1) 



5 x{2z^ + zr — r^) 



167r 



j.2p 



sin (/)(2 cos^ 



cos 61 - 1) 



5 y{2z^ + zr — r^) 
167r r'^p 



167r 
5 

47r 



sin 0(cos 6l + l)(2cos2</>- 1) 



5 (z + r)(.T^ — y^) 



j,2p 



cos ^ sin sin 6* (cos ^? + 1) = W ^ 

V 47r r-^p 



(B6a) 
(B6b) 
(B6c) 
(B6d) 
(B6e) 
(B6f) 
(B6g) 
(B6h) 
(B6i) 
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